%===================================================================================
\section{Adjoint sensitivity analysis example problems}\label{s:adj_examples}
%===================================================================================

The next three sections describe in detail a serial example
(\id{cvsRoberts\_ASAi\_dns}) and two parallel examples
(\id{cvsAdvDiff\_ASAp\_non\_p} and \id{cvsAtmDisp\_ASAi\_kry\_bbd\_p})
that perform adjoint sensitivity analysis.  For details on the other
examples, the reader is directed to the comments in their source files.

%------------------------------------------------------------------------
\subsection{A serial dense example: cvsRoberts\_ASAi\_dns}\label{ss:cvsRoberts_ASAi_dns}

As a first example of using {\cvodes} for adjoint sensitivity analysis,
we examine the chemical kinetics problem (from \id{cvsRoberts\_FSA\_dns}) 
\begin{equation}\label{e:cvsRoberts_ASAi_dns_ODE}
  \begin{split}
    &{\dot y}_1 = -p_1 y_1 + p_2 y_2 y_3   \\
    &{\dot y}_2 =  p_1 y_1 - p_2 y_2 y_3 - p_3 y_2^2 \\
    &{\dot y}_3 =  p_3 y_2^2 \\
    &y(t_0) = y_0 \, ,
  \end{split}
\end{equation}
for which we want to compute the gradient with respect to $p$ of 
\begin{equation}\label{e:cvsRoberts_ASAi_dns_G}
  G(p) = \int_{t_0}^{T}  y_3  dt ,
\end{equation}
without having to compute the solution sensitivities ${dy}/{dp}$.
Following the derivation in \ugref{ss:adj_sensi}, and taking into account
the fact that the initial values of (\ref{e:cvsRoberts_ASAi_dns_ODE}) do not depend on 
the parameters $p$, by (\ref{e:dGdp}) this gradient is simply
\begin{equation}\label{e:cvsRoberts_ASAi_dns_dGdp}
\frac{dG}{dp} = \int_{t_0}^{T} 
\left( g_p + \lambda^T f_p \right) dt \, ,
\end{equation}
where $g(t,y,p) = y_3$, $f$ is the vector-valued function 
defining the right-hand side of (\ref{e:cvsRoberts_ASAi_dns_ODE}), and $\lambda$ is 
the solution of the adjoint problem (\ref{e:adj_eqns}),
\begin{equation}\label{e:cvsRoberts_ASAi_dns_ADJ}
  \begin{split}
    &{\dot\lambda} = - (f_y)^T  \lambda - (g_y)^T \\
    &\lambda(T) = 0 \, .
  \end{split}
\end{equation}

In order to avoid saving intermediate $\lambda$ values just for the
evaluation of the integral in (\ref{e:cvsRoberts_ASAi_dns_dGdp}), we extend the
backward problem with the following $N_p$ quadrature equations
\begin{equation}\label{e:cvsRoberts_ASAi_dns_XI}
  \begin{split}
    &{\dot\xi} = g_p^T + f_p^T \lambda \\
    &\xi (T) = 0 \, ,
  \end{split}
\end{equation}
which yield $\xi(t_0) = - \int_{t_0}^{T} ( g_p^T + f_p^T \lambda) dt$
and thus ${dG}/{dp} = -\xi^T(t_0)$.
Similarly, the value of $G$ in (\ref{e:cvsRoberts_ASAi_dns_G}) can be obtained as
$G = - \zeta(t_0)$, where $\zeta$ is solution of the following quadrature
equation:
\begin{equation}\label{e:cvsRoberts_ASAi_dns_ZETA}
  \begin{split}
    &{\dot\zeta} = g = y_3 \\
    &\zeta(T) = 0 \, .
  \end{split}
\end{equation}

The main program and the user-defined routines are described below, 
with emphasis on the aspects particular to adjoint sensitivity calculations.

The calling program includes the {\cvodes} header files \id{cvodes.h} 
for {\cvodes} definitions and interface function prototypes,
the header file \id{nvector\_serial.h} for the definition of the serial implementation
of the {\nvector} module, {\nvecs}, the header files
\id{sunmatrix\_dense.h} and \id{sunlinsol\_dense.h} for the dense
{\sunmatrix} and {\sunlinsol} modules, the header file
\id{sundials\_types.h} for the definition of \id{sunrealtype} and \id{sunindextype},
and the file \id{sundials\_math.h} for the
definition of the \id{SUNRabs} macro.
%% 
This program also includes two user-defined accessor macros, \id{Ith} and \id{IJth},
that are useful in writing the problem functions in a form closely matching their
mathematical description, i.e. with components numbered from 1 instead of from 0. 
Following that, the program defines problem-specific constants and a user-defined 
data structure, which will be used to pass the values of the parameters $p$ to various
user routines. The constant \id{STEPS} defines the number of integration steps
between two consecutive checkpoints.
The program prologue ends with the prototypes of four user-supplied functions that are
called by {\cvodes}. The first two provide the right-hand side and dense Jacobian
for the forward problem, and the last two provide the right-hand side and dense
Jacobian for the backward problem.

The \id{main} function begins with type declarations and continues with the
allocation and initialization of the user data structure, which contains the values 
of the parameters $p$. Next, it allocates and
initializes \id{y} with the initial conditions for the forward problem, allocates and 
initializes \id{q} for the quadrature used in computing the value $G$, and finally
sets the scalar relative tolerance \id{reltolQ} and vector absolute tolerance
\id{abstolQ} for the quadrature variables.
No tolerances for the state variables are defined since \id{cvsRoberts\_ASAi\_dns}
uses its own function to compute the error weights for WRMS norm
estimates of state solution vectors.

The call to \id{CVodeCreate} creates the main integrator memory block for the 
forward integration and specifies the \id{CV\_BDF} integration method. 
The call to \id{CVodeInit} initializes the forward integration by specifying the 
initial conditions.
The call to \id{CVodeWFtolerances} specifies a function that computes error weights.
The next call specifies the optional user data pointer \id{data}.
The linear solver is selected to be {\sunlinsoldense} through calls to
create the template Jacobian matrix and dense linear solver objects
(\id{SUNDenseMatrix} and \id{SUNLinSol\_Dense}), and to attach these
to the {\cvodes} integrator via the call to
\id{CVodeSetLinearSolver}. The user-provided Jacobian routine \id{Jac}
is specified through a call to \id{CVodeSetJacFn}. 

The next code block initializes quadrature computations in the forward phase, by
allocating {\cvodes} memory for quadrature integration (the call to
\id{CVodeQuadInit} specifies the right-hand side \id{fQ} of the quadrature
equation and the initial values of the quadrature variable),
setting the integration tolerances for the quadrature variables, and finally
including the quadrature variable in the error test.

Allocation for the memory block of the combined forward-backward
problem is accomplished through the call to \id{CVadjInit} which
specifies \id{STEPS = 150}, the number of steps between two
checkpoints, and specifies cubic Hermite interpolation.

The call to \id{CVodeF} requests the solution of the forward problem to \id{TOUT}.
If successful, at the end of the integration, \id{CVodeF} will return the number
of saved checkpoints in the argument \id{ncheck} (optionally, a list of the
checkpoints can be obtained by calling \id{CVodeGetAdjCheckPointsInfo} and the
checkpoint information printed).

The next segment of code deals with the setup of the backward problem. 
First, a serial vector \id{yB} of length \id{NEQ} is allocated and initialized with the
value of $\lambda (= 0.0)$ at the final time (\id{TB1} = 4.0E7).  A second
serial vector \id{qB} of dimension \id{NP} is created and initialized to $0.0$.
This vector corresponds to the quadrature variables $\xi$ whose values at $t_0$
will be the components of the desired gradient of $\partial G / \partial p$
(after a sign change).
%%
Following that, the program sets the relative and absolute tolerances
for the backward integration.

The {\cvodes} memory for the backward integration is created and allocated
by the calls to the interface routines \id{CVodeCreateB} and \id{CVodeInitB} which 
specify the \id{CV\_BDF} integration method, among other things.  The
dense linear solver is created and initialized by calling the
\id{SUNDenseMatrix}, \id{SUNLinSol\_Dense} and
\id{CVodeSetLinearSolverB} routines, and specifying a non-\id{NULL}
Jacobian routine \id{JacB} and user data \id{data}.

The tolerances for the integration of quadrature variables, \id{reltolB} and
\id{abstolQB}, are specified through \id{CVodeQuadSStolerancesB}.
The call to \id{CVodeSetQuadErrConB} indicates that $\xi$ should be included
in the error test.
Quadrature computation is initialized by calling \id{CVodeQuadInitB}
which specifies the right-hand side of the quadrature equations as \id{fQB}.

The actual solution of the backward problem is accomplished through
two calls to \id{CVodeB} --- one for intermediate output at $t = 40$,
and one for the final time \id{T0} $= 0$.  At each point, the backward
solution \id{yB} ($ = \lambda$) is obtained with a call to \id{CVodeGetB} and the
forward solution with a call to \id{CVodeGetAdjY}.  The values of the
quadrature variables $\xi$ at time \id{T0} are loaded in \id{qB} by calling
the extraction routine \id{CVodeGetQuadB}. The negative of \id{qB}
gives the gradient $\partial G / \partial p$.

The main program then carries out a second backward problem.
It calls to \id{CVodeReInitB} and \id{CVodeQuadReInitB} to re-initialize the 
backward memory block for a new adjoint computation with a different final 
time (\id{TB2} $= 50$).  This is followed by two calls to \id{CVodeB},
one for intermediate output at $t = 40$ and one for the final values at
$t = 0$.  Finally, the gradient $\partial G / \partial p$ of the second
function $G$ is printed.

The main program ends by freeing previously allocated memory by calling 
\id{CVodeFree} (for the {\cvodes} memory for the forward problem), \id{CVadjFree} 
(for the memory allocated for the combined problem), and \id{N\_VFree\_Serial} 
(for the various vectors).

The user-supplied functions \id{f} and \id{Jac} for the right-hand side and
Jacobian of the forward problem are straightforward expressions of its 
mathematical formulation (\ref{e:cvsRoberts_ASAi_dns_ODE}). 
The function \id{ewt} is the same as the one for \id{cvRoberts\_dns\_uw.c}.
The function \id{fQ} implements
(\ref{e:cvsRoberts_ASAi_dns_ZETA}), while \id{fB}, \id{JacB}, and \id{fQB} are mere
translations of the backward problem (\ref{e:cvsRoberts_ASAi_dns_ADJ}) and
(\ref{e:cvsRoberts_ASAi_dns_XI}).

The output generated by \id{cvsRoberts\_ASAi\_dns} is shown below.

\includeOutput{cvsRoberts\_ASAi\_dns}{../../examples/cvodes/serial/cvsRoberts_ASAi_dns.out}

%--------------------------------------------------------------------------

\subsection{A parallel nonstiff example: cvsAdvDiff\_ASAp\_non\_p}
\label{ss:cvsAdvDiff_ASAp_non_p}

As an example of using the {\cvodes} adjoint sensitivity module with the
parallel vector module {\nvecp}, we describe a sample program that solves
the following problem:  Consider the 1-D advection-diffusion equation
\begin{equation}\label{e:cvsAdvDiff_ASAp_non_p:orig_pde}
  \begin{split}
    & \frac{\partial u}{\partial t} = p_1 \frac{\partial^2 u}{\partial x^2} 
    + p_2 \frac{\partial u}{\partial x} \\
    & 0 = x_0 \le x \le x_1 = 2 \\
    & 0 = t_0 \le t \le t_f = 2.5 \, ,
  \end{split}
\end{equation}
with boundary conditions $u(t,x_0) = u(t,x_1) = 0 ,\, \forall t$,
and initial condition $u(t_0 , x) = u_0(x) = x(2-x)e^{2x}$. Also
consider the function
\begin{equation*}
  g(t) = \int_{x_0}^{x_1} u(t,x) dx \, .
\end{equation*}
We wish to find, through adjoint sensitivity analysis, the gradient of
$g(t_f)$ with respect to $p = [p_1 ; p_2]$ and the perturbation in $g(t_f)$
due to a perturbation $\delta u_0$ in $u_0$.

The approach we take in the program \id{cvsAdvDiff\_ASAp\_non\_p} is to first derive an 
adjoint PDE which is then discretized in space and integrated backwards
in time to yield the desired sensitivities. A straightforward extension 
to PDEs of the derivation given in \ugref{ss:adj_sensi} gives
\begin{equation}\label{e:cvsAdvDiff_ASAp_non_p:dgdp}
  \frac{dg}{dp} (t_f) = \int_{t_0}^{t_f} dt 
  \int_{x_0}^{x_1} dx \mu \cdot 
  \left[
    \frac{\partial^2 u}{\partial x^2} ;
    \frac{\partial u}{\partial x}
  \right ]
\end{equation}
and
\begin{equation}\label{e:cvsAdvDiff_ASAp_non_p:delg}
  \delta g |_{t_f} = \int_{x_0}^{x_1} \mu(t_0,x) \delta u_0(x) dx \, , 
\end{equation}
where $\mu$ is the solution of the adjoint PDE
\begin{equation}\label{e:cvsAdvDiff_ASAp_non_p:adj_pde}
  \begin{split}
    & \frac{\partial \mu}{\partial t} + p_1 \frac{\partial^2 \mu}{\partial x^2} 
    - p_2 \frac{\partial \mu}{\partial x} = 0 \\
    & \mu(t_f , x) = 1 \\
    & \mu(t , x_0) = \mu( t , x_1 ) = 0 \, .
  \end{split}
\end{equation}
Both the forward problem (\ref{e:cvsAdvDiff_ASAp_non_p:orig_pde}) and the backward
problem (\ref{e:cvsAdvDiff_ASAp_non_p:adj_pde}) are discretized on a uniform spatial
grid of size $M_x + 2$ with central differencing and with boundary values eliminated,
leaving ODE systems of size $N = M_x$ each.  As always, we deal with the time
quadratures in (\ref{e:cvsAdvDiff_ASAp_non_p:dgdp}) by introducing the additional
equations
\begin{equation}\label{e:cvsAdvDiff_ASAp_non_p:quad}
  \begin{split}
    &{\dot\xi}_1 = \int_{x_0}^{x_1} dx \mu \frac{\partial^2 u}{\partial x^2} \, , \quad
    \xi_1(t_f) = 0 \, , \\
    &{\dot\xi}_2 = \int_{x_0}^{x_1} dx \mu \frac{\partial u}{\partial x} \, , \quad
    \xi_2(t_f) = 0 \, ,
  \end{split}
\end{equation}
yielding
\begin{equation*}
  \frac{dg}{dp} (t_f) = - \left[ \xi_1(t_0) ; \xi_2(t_0) \right ]
\end{equation*}
The space integrals in (\ref{e:cvsAdvDiff_ASAp_non_p:delg}) and
(\ref{e:cvsAdvDiff_ASAp_non_p:quad}) are
evaluated numerically, on the given spatial mesh, using the trapezoidal rule.

Note that $\mu(t_0 , x^*)$ is nothing but the perturbation in $g(t_f)$
due to a $\delta$-function perturbation $\delta u_0(x) = \delta(x-x^*)$ in the
initial conditions.  Therefore, $\mu(t_0,x)$ completely describes $\delta g(t_f)$
for any perturbation $\delta u_0$.

Both the forward and the backward problems are solved with the option for nonstiff systems,
i.e. using the Adams method with fixed-point iteration for the solution of
the nonlinear systems. The overall structure of the \id{main} function is very
similar to that of the code \id{cvsRoberts\_ASAi\_dns} discussed previously with 
differences arising from the use of the parallel {\nvector} module. Unlike 
\id{cvsRoberts\_ASAi\_dns}, the example \id{cvsAdvDiff\_ASAp\_non\_p} illustrates
computation of the additional quadrature variables by appending \id{NP} equations
to the adjoint system.  This approach can be a better alternative to using special
treatment of the quadrature equations when their number is too small for parallel
treatment.

Besides the parallelism implemented by {\cvodes} at the {\nvector} level,
this example uses {\mpi} calls to parallelize the calculations of the
right-hand side routines \id{f} and \id{fB} and of the spatial integrals involved.
The forward problem has size \id{NEQ = MX}, while the backward problem has
size \id{NB = NEQ + NP}, where \id{NP = 2} is the number of quadrature equations
in (\ref{e:cvsAdvDiff_ASAp_non_p:quad}).
The use of the total number of available processes on two problems of different 
sizes deserves some comments, as this is typical in adjoint sensitivity 
analysis. Out of the total number of available processes, namely \id{nprocs},
the first \id{npes = nprocs - 1} processes are dedicated to the integration of
the ODEs arising from the semi-discretization of the PDEs 
(\ref{e:cvsAdvDiff_ASAp_non_p:orig_pde}) and (\ref{e:cvsAdvDiff_ASAp_non_p:adj_pde}),
and receive the same load on both the forward and backward integration phases. 
The last process is reserved for the integration of the quadrature equations 
(\ref{e:cvsAdvDiff_ASAp_non_p:quad}), and is therefore inactive during the forward 
phase. Of course, for problems involving a much larger number of quadrature equations,
more than one process could be reserved for their integration. 
An alternative would be to redistribute the \id{NB} backward problem variables 
over all available processes, without any relationship to the load distribution 
of the forward phase. However, the approach taken in \id{cvsAdvDiff\_ASAp\_non\_p}
has the advantage that the communication strategy adopted for the forward problem 
can be directly transferred to communication among the first \id{npes}
processes during the backward integration phase. 

We must also emphasize that, although inactive during the forward integration phase, 
the last process {\em must} participate in that phase with a 
{\em zero local array length}. 
This is because, during the backward integration phase, this process must
have its own local copy of variables (such as \id{cvadj\_mem}) that were set
only during the forward phase.

Using \id{MX} $=40$ on 4 processes, the gradient of $g(t_f)$ with respect to 
the two problem parameters is obtained as $dg/dp(t_f) = [ -1.13856; -1.01023]$.
The gradient of $g(t_f)$ with respect to the initial conditions is shown in
Fig.~\ref{f:cvsAdvDiff_ASAp_non_p}. The gradient is plotted superimposed over the
initial conditions.
%%
\begin{figure}
  {\centerline{\includegraphics[width=0.75\textwidth]{cvsadjnonx_p}}}
  \caption{Results for the \id{cvsAdvDiff\_ASAp\_non\_p} example problem.
    The gradient of $g(t_f)$ with respect to the initial conditions $u_0$ 
    is shown superimposed over the values $u_0$.}
  \label{f:cvsAdvDiff_ASAp_non_p}
\end{figure}
%%
Sample output generated by \id{cvsAdvDiff\_ASAp\_non\_p}, for \id{MX} $=20$, is
shown below.

\includeOutput{cvsAdvDiff\_ASAp\_non\_p}{../../examples/cvodes/parallel/cvsAdvDiff_ASAp_non_p.out}


%--------------------------------------------------------------------------

\subsection{A parallel example using CVBBDPRE: cvsAtmDisp\_ASAi\_kry\_bbd\_p}
\label{ss:cvsAtmDisp_ASAi_kry_bbd_p}

As a more elaborate example of a parallel adjoint sensitivity calculation,
we describe next the program \id{cvsAtmDisp\_ASAi\_kry\_bbd\_p} provided
with {\cvodes}. This example models an atmospheric 
release with an advection-diffusion PDE in 2-D or 3-D and computes the gradient 
with respect to source parameters of the space-time average of the squared norm
of the concentration.
Given a known velocity field $v(t,x)$ and source function $S$, the  transport
equation for the concentration $c(t,x)$ in a domain $\Omega$ is given by
\begin{equation}\label{e:cvsAtmDisp_ASAi_kry_bbd_p_PDE}
  \begin{split}
    \frac{\partial c}{\partial t} - k \nabla^2 c + v \cdot \nabla c + S = 0 \, , 
    &\text{ in } (0,T) \times \Omega \\
    \frac{\partial c}{\partial n} = g \, ,
    &\text{ on } (0,T) \times \partial\Omega \\
    c = c_0(x) \, ,
    &\text{ in } \Omega \text{ at } t = 0 \, ,
  \end{split}
\end{equation}
where $\Omega$ is a box in ${\mathbb{R}}^2$ or ${\mathbb{R}}^3$ and $n$ is the 
normal to the boundary of $\Omega$.
We assume homogeneous boundary conditions ($g = 0$) and a zero initial
concentration everywhere in $\Omega$ ($c_0(x) = 0$). The wind field has only a
nonzero component in the $x$ direction given by a Poiseuille profile along the 
direction $y$.

Using adjoint sensitivity analysis, the gradient of
\begin{equation}\label{e:cvsAtmDisp_ASAi_kry_bbd_p_G}
  G(p) = \frac{1}{2} \int_0^T \int_\Omega \| c(t,x) \|^2 \, d\Omega \, dt
\end{equation}
is obtained as
\begin{equation}\label{e:cvsAtmDisp_ASAi_kry_bbd_p_dGdp}
  \frac{dG}{dp_i} = \int_t \int_\Omega \lambda(t,x) \delta(x-x_i) \, d\Omega \, dt
  = \int_t \lambda(t,x_i) \, dt \, ,
\end{equation}
where $x_i$ is the location of the source of intensity $S(x_i)=p_i$, and $\lambda$
is solution of the adjoint PDE
\begin{equation}\label{e:cvsAtmDisp_ASAi_kry_bbd_p_ADJ}
  \begin{split}
    - \frac{\partial\lambda}{\partial t} - k \nabla^2\lambda - v \cdot \lambda = c(t,x)  \, ,
    &\text{ in } (T,0) \times \Omega \\
    (k \nabla\lambda + v \lambda) \cdot n = 0 \, ,
    &\text{ on } (0,T) \times \partial\Omega \\
    \lambda = 0 \, ,
    &\text{ in } \Omega \text{ at } t = T \, .
  \end{split}
\end{equation}
%%
The PDE (\ref{e:cvsAtmDisp_ASAi_kry_bbd_p_PDE}) is semi-discretized in space with 
central finite differences, with the boundary conditions explicitly taken into account 
by using layers of ghost cells in every direction. If the direction $x^i$ of $\Omega$
is discretized into $m_i$ intervals, this leads to a system of ODEs of dimension 
$N = \prod_1^d (m_i+1)$, with $d=2$, or $d=3$.
The source term $S$ is parameterized as a piecewise constant function and yielding
$N$ parameters in the problem. The nominal values of the source parameters correspond
to two Gaussian sources.

The source code as supplied runs the 2-D problem.  To obtain the 3-D version,
add a line \id{\#define~USE3D} at the top of \id{main}.

The adjoint PDE (\ref{e:cvsAtmDisp_ASAi_kry_bbd_p_ADJ}) is discretized to a system of
ODEs in a similar fashion.  The space integrals in (\ref{e:cvsAtmDisp_ASAi_kry_bbd_p_G})
and (\ref{e:cvsAtmDisp_ASAi_kry_bbd_p_dGdp}) are simply approximated by their
Riemann sums, while the time integrals are resolved by appending pure quadrature
equations to the systems of ODEs.

We use BDF with the {\sunlinsolspgmr} linear solver module and the
{\cvbbdpre} preconditioner for both the forward and the backward
integration phases. The value of $G$ is computed on the forward phase
as a quadrature, while the components of the gradient $dG/dp$ are
computed as quadratures during the backward integration phase.  All
quadrature variables are included in the corresponding error tests.

Communication between processes for the evaluation of the ODE right-hand sides involves
passing the solution on the local boundaries (lines in 2-D, surfaces in 3-D) to 
the 4 (6 in 3-D) neighboring processes. This is implemented in the function 
\id{f\_comm}, called in \id{f} and \id{fB} before evaluation of the local residual 
components. Since there is no additional communication required for the {\cvbbdpre}
preconditioner, a \id{NULL} pointer is passed for \id{gloc} and \id{glocB} in the
calls to \id{CVBBDPrecInit} and \id{CVBBDPrecInitB}, respectively.

For the sake of clarity, the \id{cvsAtmDisp\_ASAi\_kry\_bbd\_p} example does not use
the most memory-efficient implementation possible, as the local segment of the 
solution vectors (\id{y} on the forward phase and \id{yB} on the backward phase)
and the data received from neighboring processes is loaded into a temporary 
array \id{y\_ext} which is then used exclusively in computing the local components
of the right-hand sides.

Note that if \id{cvsAtmDisp\_ASAi\_kry\_bbd\_p} is given any command line argument,
it will generate a series of MATLAB files which can be used to visualize the solution.
The results of a 2-D simulation and adjoint sensitivity analysis with
\id{cvsAtmDisp\_ASAi\_kry\_bbd\_p} on a $80 \times 80$ grid and $2 \times 4 = 8$
processes are shown in Fig.~\ref{f:cvsAtmDisp_ASAi_kry_bbd_p2D}.
Results in 3-D\footnote{The name of the executable for the 3-D version is
\id{cvsAtmDisp\_ASAi\_kry\_bbd\_p3D}.},  on a $80 \times 80 \times 40$ grid and 
$2 \times 4 \times 2= 16$ processes are shown in 
Figs.~\ref{f:cvsAtmDisp_ASAi_kry_bbd_p3D_a} and \ref{f:cvsAtmDisp_ASAi_kry_bbd_p3D_b}.
%%
\begin{figure}
  {\centerline{\includegraphics[width=\textwidth]{cvsadjkryx_p2D}}}
  \caption{Results for the \id{cvsAtmDisp\_ASAi\_kry\_bbd\_p} example problem in 2D. 
    The gradient with respect to the source parameters is pictured on the left. 
    On the right, the gradient was color-coded and superimposed over the nominal value 
    of the source parameters.}
  \label{f:cvsAtmDisp_ASAi_kry_bbd_p2D}
\end{figure}
%%
%%
\begin{figure}
  {\centerline{\includegraphics[width=0.75\textwidth]{cvsadjkryx_p3Dcf.png}}}
  \caption{Results for the \id{cvsAtmDisp\_ASAi\_kry\_bbd\_p} example problem in 3D.
  Nominal values of the source parameters.}\label{f:cvsAtmDisp_ASAi_kry_bbd_p3D_a}
\end{figure}
%%
\begin{figure}
  {\centerline{\includegraphics[width=0.75\textwidth]{cvsadjkryx_p3Dgrad.png}}}
  \caption{Results for the \id{cvsAtmDisp\_ASAi\_kry\_bbd\_p} example problem in 3D.
  Two isosurfaces of the gradient with respect to the source parameters. They correspond
  to values of $0.25$ (green) and $0.4$ (blue).}\label{f:cvsAtmDisp_ASAi_kry_bbd_p3D_b}  
\end{figure}
%%
A sample output generated by \id{cvsAtmDisp\_ASAi\_kry\_bbd\_p} for a 2D calculation
is shown below.

\includeOutput{cvsAtmDisp\_ASAi\_kry\_bbd\_p}{../../examples/cvodes/parallel/cvsAtmDisp_ASAi_kry_bbd_p.out}

